Code
library(here)
library(flextable)
library(plotly)
source(here("02_scripts/03a_fun_calib_td.R"))
source(here("02_scripts/06_interventions_models.R"))
theme_set(theme_minimal())05_cba_results
Scenario 1: Accounts for the effect of Fentanyl contamination starting from 2018 Scenario 2: Scenario 2: Doesn’t account for Fentanyl contamination and has the same probability from Illicit use to OD and from OD to OD death
Code
trace_tbl <- as_tibble(mod_basecase$m_M) %>%
mutate(cycle_num = 0:n_cycles) %>%
pivot_longer(1:n_states, names_to = "state", values_to = "count") %>%
mutate(mod = "No Interventions") %>%
bind_rows(., as_tibble(mod_nalx_1$m_M) %>%
mutate(cycle_num = 0:n_cycles) %>%
pivot_longer(1:n_states, names_to = "state", values_to = "count") %>%
mutate(mod = "Naloxone")) %>%
bind_rows(., as_tibble(mod_ss_1$m_M) %>%
mutate(cycle_num = 0:n_cycles) %>%
pivot_longer(1:n_states, names_to = "state", values_to = "count") %>%
mutate(mod = "Safer Supply")) %>%
bind_rows(., as_tibble(mod_pres_guid_1$m_M) %>%
mutate(cycle_num = 0:n_cycles) %>%
pivot_longer(1:n_states, names_to = "state", values_to = "count") %>%
mutate(mod = "Prescription Guidelines")) %>%
bind_rows(., as_tibble(mod_all_int_1$m_M) %>%
mutate(cycle_num = 0:n_cycles) %>%
pivot_longer(1:n_states, names_to = "state", values_to = "count") %>%
mutate(mod = "All Interventions"))
trace_tbl$mod <- factor(trace_tbl$mod,
levels = c("No Interventions", "Naloxone", "Safer Supply", "Prescription Guidelines", "All Interventions"),
labels = c("No Interventions", "Naloxone", "Safer Supply", "Prescription Guidelines", "All Interventions"))
v_extra_sevbi <- c(mod_basecase$extra_sev_bi, mod_nalx_1$extra_sev_bi,
mod_ss_1$extra_sev_bi, mod_pres_guid_1$extra_sev_bi,
mod_all_int_1$extra_sev_bi)
v_extra_modbi <- c(mod_basecase$extra_mod_bi, mod_nalx_1$extra_mod_bi,
mod_ss_1$extra_mod_bi, mod_pres_guid_1$extra_mod_bi,
mod_all_int_1$extra_mod_bi)
v_extra_od <- c(mod_basecase$extra_od_deaths, mod_nalx_1$extra_od_deaths,
mod_ss_1$extra_od_deaths, mod_pres_guid_1$extra_od_deaths,
mod_all_int_1$extra_od_deaths)
mod_names <- c("No Interventions", "Naloxone", "Safer Supply", "Prescription Guidelines", "All Interventions")
mod_names <- factor(mod_names, levels = mod_names, labels = mod_names)
for (j in 1:length(mod_names)){
trace_tbl$count[trace_tbl$state == "BO_OD_DEATH" &
trace_tbl$cycle_num == 180 &
trace_tbl$mod == mod_names[j]] <- trace_tbl$count[trace_tbl$state == "BO_OD_DEATH" &
trace_tbl$cycle_num == 180 &
trace_tbl$mod == mod_names[j]] + v_extra_od[j]
trace_tbl$count[trace_tbl$state == "BO_MOD_BI" &
trace_tbl$cycle_num == 180 &
trace_tbl$mod == mod_names[j]] <- trace_tbl$count[trace_tbl$state == "BO_MOD_BI" &
trace_tbl$cycle_num == 180 &
trace_tbl$mod == mod_names[j]] + v_extra_modbi[j]
trace_tbl$count[trace_tbl$state == "BO_SEVERE_BI" &
trace_tbl$cycle_num == 180 &
trace_tbl$mod == mod_names[j]] <- trace_tbl$count[trace_tbl$state == "BO_SEVERE_BI" &
trace_tbl$cycle_num == 180 &
trace_tbl$mod == mod_names[j]] + v_extra_sevbi[j]
}
trace_tbl$state <- factor(trace_tbl$state,
levels = v_state_names,
labels = v_state_names)Code
trace_tbl_2 <- as_tibble(mod_basecase_2$m_M) %>%
mutate(cycle_num = 0:n_cycles) %>%
pivot_longer(1:n_states, names_to = "state", values_to = "count") %>%
mutate(mod = "No Interventions") %>%
bind_rows(., as_tibble(mod_nalx_2$m_M) %>%
mutate(cycle_num = 0:n_cycles) %>%
pivot_longer(1:n_states, names_to = "state", values_to = "count") %>%
mutate(mod = "Naloxone")) %>%
bind_rows(., as_tibble(mod_ss_2$m_M) %>%
mutate(cycle_num = 0:n_cycles) %>%
pivot_longer(1:n_states, names_to = "state", values_to = "count") %>%
mutate(mod = "Safer Supply")) %>%
bind_rows(., as_tibble(mod_pres_guid_2$m_M) %>%
mutate(cycle_num = 0:n_cycles) %>%
pivot_longer(1:n_states, names_to = "state", values_to = "count") %>%
mutate(mod = "Prescription Guidelines")) %>%
bind_rows(., as_tibble(mod_all_int_2$m_M) %>%
mutate(cycle_num = 0:n_cycles) %>%
pivot_longer(1:n_states, names_to = "state", values_to = "count") %>%
mutate(mod = "All Interventions"))
trace_tbl_2$mod <- factor(trace_tbl_2$mod,
levels = c("No Interventions", "Naloxone", "Safer Supply", "Prescription Guidelines", "All Interventions"),
labels = c("No Interventions", "Naloxone", "Safer Supply", "Prescription Guidelines", "All Interventions"))
v_extra_sevbi_2 <- c(mod_basecase_2$extra_sev_bi, mod_nalx_2$extra_sev_bi,
mod_ss_2$extra_sev_bi, mod_pres_guid_2$extra_sev_bi,
mod_all_int_2$extra_sev_bi)
v_extra_modbi_2 <- c(mod_basecase_2$extra_mod_bi, mod_nalx_2$extra_mod_bi,
mod_ss_2$extra_mod_bi, mod_pres_guid_2$extra_mod_bi,
mod_all_int_2$extra_mod_bi)
v_extra_od_2 <- c(mod_basecase_2$extra_od_deaths, mod_nalx_2$extra_od_deaths,
mod_ss_2$extra_od_deaths, mod_pres_guid_2$extra_od_deaths,
mod_all_int_2$extra_od_deaths)
for (j in 1:length(mod_names)){
trace_tbl_2$count[trace_tbl_2$state == "BO_OD_DEATH" &
trace_tbl_2$cycle_num == 180 &
trace_tbl_2$mod == mod_names[j]] <- trace_tbl_2$count[trace_tbl_2$state == "BO_OD_DEATH" &
trace_tbl_2$cycle_num == 180 &
trace_tbl_2$mod == mod_names[j]] + v_extra_od_2[j]
trace_tbl_2$count[trace_tbl_2$state == "BO_MOD_BI" &
trace_tbl_2$cycle_num == 180 &
trace_tbl_2$mod == mod_names[j]] <- trace_tbl_2$count[trace_tbl_2$state == "BO_MOD_BI" &
trace_tbl_2$cycle_num == 180 &
trace_tbl_2$mod == mod_names[j]] + v_extra_modbi_2[j]
trace_tbl_2$count[trace_tbl_2$state == "BO_SEVERE_BI" &
trace_tbl_2$cycle_num == 180 &
trace_tbl_2$mod == mod_names[j]] <- trace_tbl_2$count[trace_tbl_2$state == "BO_SEVERE_BI" &
trace_tbl_2$cycle_num == 180 &
trace_tbl_2$mod == mod_names[j]] + v_extra_sevbi_2[j]
}
trace_tbl_2$state <- factor(trace_tbl_2$state,
levels = v_state_names,
labels = v_state_names)
trace_tbl <- trace_tbl %>%
mutate(scenario = "Scenario 1") %>%
bind_rows(., trace_tbl_2 %>%
mutate(scenario = "Scenario 2"))1 Cumulative OD Deaths
1.0.0.1 Cumulative OD Deaths - Scenario 1
Code
p1 <- ggplot(trace_tbl %>%
filter(state == "BO_OD_DEATH") %>%
filter(cycle_num != 0) %>%
filter((cycle_num + 1) %% 12 == 1) %>%
filter(scenario == "Scenario 1") %>%
mutate(year = rep(c(2015:2029), length(mod_names))),
aes(x = year, y = count, colour = mod)) +
geom_line() +
geom_point(size = 1) +
ylab("Cumulative OD Deaths") + xlab("Year") +
labs(colour = "Intervention")
ggplotly(p1)1.0.0.2 Cumulative OD Deaths - Scenario 2
Code
p2 <- ggplot(trace_tbl %>%
filter(state == "BO_OD_DEATH") %>%
filter(cycle_num != 0) %>%
filter((cycle_num + 1) %% 12 == 1) %>%
filter(scenario == "Scenario 2") %>%
mutate(year = rep(c(2015:2029), length(mod_names))),
aes(x = year, y = count, colour = mod)) +
geom_line() +
geom_point(size = 1) +
ylab("Cumulative OD Deaths") + xlab("Year") +
labs(colour = "Intervention")
ggplotly(p2)1.0.0.3 Both scenarios
Code
plotly::ggplotly(ggplot(trace_tbl %>%
filter(state == "BO_OD_DEATH") %>%
filter(cycle_num != 0) %>%
filter((cycle_num + 1) %% 12 == 1) %>%
mutate(year = rep(c(2015:2029), length(mod_names)*2)),
aes(x = year, y = count, colour = mod)) +
geom_line() +
geom_point(size = 1) +
ylab("Cumulative OD Deaths") + xlab("Year") +
labs(colour = "Intervention") +
facet_wrap(~scenario) +
labs(caption = "Scenario 1: Accounts for the effect of Fentanyl contamination starting from 2018 \n
Scenario 2: Doesn't account for Fentanyl contamination and has the same probability from \n Illicit use to OD and from OD to OD death"))2 New OD Deaths per year
Code
inci_od_deaths_basecase <- inci_od_deaths_nalox <- inci_od_deaths_ss <- inci_od_deaths_pg <- inci_od_deaths_all <- inci_od_deaths_basecase_2 <- inci_od_deaths_nalox_2 <- inci_od_deaths_ss_2 <- inci_od_deaths_pg_2 <- inci_od_deaths_all_2 <- rep(NA, 14)
for (i in 1:14){
yr <- 2015 + i
inci_od_deaths_basecase[i] <- mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
inci_od_deaths_basecase_2[i] <- mod_basecase_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_basecase_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
inci_od_deaths_nalox[i] <- mod_nalx_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_nalx_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
inci_od_deaths_nalox_2[i] <- mod_nalx_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_nalx_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
inci_od_deaths_ss[i] <- mod_ss_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_ss_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
inci_od_deaths_ss_2[i] <- mod_ss_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_ss_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
inci_od_deaths_pg[i] <- mod_pres_guid_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_pres_guid_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
inci_od_deaths_pg_2[i] <- mod_pres_guid_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_pres_guid_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
inci_od_deaths_all[i] <- mod_all_int_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_all_int_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
inci_od_deaths_all_2[i] <- mod_all_int_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_all_int_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
}
inc_od_death_tbl <- data.frame(intervention = c(rep(mod_names, each = 14)),
year = c(rep(2016:2029, 5)),
inci_od_deaths = c(inci_od_deaths_basecase,
inci_od_deaths_nalox,
inci_od_deaths_ss,
inci_od_deaths_pg,
inci_od_deaths_all),
scenario = rep("Scenario 1", 14*5)) %>%
bind_rows(., data.frame(intervention = c(rep(mod_names, each = 14)),
year = c(rep(2016:2029, 5)),
inci_od_deaths = c(inci_od_deaths_basecase_2,
inci_od_deaths_nalox_2,
inci_od_deaths_ss_2,
inci_od_deaths_pg_2,
inci_od_deaths_all_2),
scenario = rep("Scenario 2", 14*5)))
saveRDS(inc_od_death_tbl, file = "inc_od_death_tbl_mod.RDS")
target_od_deaths <- calib_target_tbl %>%
filter(group == "total_od_deaths") %>%
select(year, target) %>%
rename(inci_od_deaths = target) %>%
mutate(intervention = "Target")
saveRDS(target_od_deaths, file = "inc_od_death_tbl_target.RDS")
inc_od_death_tbl$intervention <- factor(inc_od_death_tbl$intervention,
levels = mod_names, labels = mod_names)2.0.0.1 Scenario 1
Code
p3 <- ggplot(inc_od_death_tbl %>% filter(scenario == "Scenario 1"),
aes(x = year, y = inci_od_deaths, color = intervention)) +
geom_line() +
geom_point(size = 1) +
geom_point(data = target_od_deaths, aes (x = year, y = inci_od_deaths)) +
ylab("New OD Deaths per year") + xlab("Year") +
labs(colour = "Intervention") +
scale_color_brewer(palette = "Dark2")
ggplotly(p3)2.0.0.2 Scenario 2
Code
p4 <- ggplot(inc_od_death_tbl %>% filter(scenario == "Scenario 2"),
aes(x = year, y = inci_od_deaths, color = intervention)) +
geom_line() +
geom_point(size = 1) +
geom_point(data = target_od_deaths, aes (x = year, y = inci_od_deaths)) +
ylab("New OD Deaths per year") + xlab("Year") +
labs(colour = "Intervention") +
scale_color_brewer(palette = "Dark2")
ggplotly(p4)2.0.0.3 Both scenarios
Code
plotly::ggplotly(ggplot(inc_od_death_tbl,
aes(x = year, y = inci_od_deaths, color = intervention)) +
geom_line() +
geom_point(size = 1) +
geom_point(data = target_od_deaths, aes (x = year, y = inci_od_deaths)) +
ylab("New OD Deaths per year") + xlab("Year") +
labs(colour = "Intervention") +
scale_color_brewer(palette = "Dark2") +
facet_wrap(~scenario))3 Cumulative Deaths
3.0.0.1 Scenario 1
Code
p5 <- ggplot(trace_tbl %>%
filter(state == "BO_DEATH") %>%
filter(cycle_num != 0) %>%
filter((cycle_num + 1) %% 12 == 1) %>%
filter(scenario == "Scenario 1") %>%
mutate(year = rep(c(2015:2029), length(mod_names))),
aes(x = year, y = count, colour = mod)) +
geom_line() +
ylab("Cumulative Deaths (excluding OD-deaths)") + xlab("Year") +
labs(colour = "Intervention")
ggplotly(p5)3.0.0.2 Scenario 2
Code
p6 <- ggplot(trace_tbl %>%
filter(state == "BO_DEATH") %>%
filter(cycle_num != 0) %>%
filter((cycle_num + 1) %% 12 == 1) %>%
filter(scenario == "Scenario 2") %>%
mutate(year = rep(c(2015:2029), length(mod_names))),
aes(x = year, y = count, colour = mod)) +
geom_line() +
ylab("Cumulative Deaths (excluding OD-deaths)") + xlab("Year") +
labs(colour = "Intervention")
ggplotly(p6)3.0.0.3 Both scenarios
Code
plotly::ggplotly(ggplot(trace_tbl %>%
filter(state == "BO_DEATH") %>%
filter(cycle_num != 0) %>%
filter((cycle_num + 1) %% 12 == 1) %>%
mutate(year = rep(c(2015:2029), length(mod_names)*2)),
aes(x = year, y = count, colour = mod)) +
geom_line() +
ylab("Cumulative Deaths (excluding OD-deaths)") + xlab("Year") +
labs(colour = "Intervention") +
facet_wrap(~scenario))4 Trace Plots
Code
set.seed(122345)
colors_pal <- Polychrome::createPalette(n_states,
c("#084C61", "#DB504A",
"#E3B505", "#4F6D7A",
"#56A3A6"))
names(colors_pal) <- NULL
plotly::ggplotly(ggplot(trace_tbl %>%
filter(grepl("BPO", state)), aes(x = cycle_num, y = count, colour = state)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
ggrepel::geom_text_repel(data = trace_tbl %>% filter(cycle_num == 180) %>%
filter(grepl("BPO", state)),
aes(x = cycle_num, y = count, label = state),
check_overlap = T, size = 3) +
scale_color_manual(values = colors_pal) +
# facet_wrap(~mod) +
facet_grid(rows = vars(mod), cols = vars(scenario),
switch = "y"))Code
plotly::ggplotly(ggplot(trace_tbl %>%
filter(state %in% c("BPO_MISUSE", "BI_ILLICIT",
"BS_OAT_INI", "BS_OAT_MAINT")), aes(x = cycle_num, y = count, colour = state)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
ggrepel::geom_text_repel(data = trace_tbl %>%
filter(state %in% c("BPO_MISUSE", "BI_ILLICIT",
"BS_OAT_INI", "BS_OAT_MAINT")) %>% filter(cycle_num == 180),
aes(x = cycle_num, y = count, label = state),
check_overlap = T, size = 3) +
scale_color_manual(values = colors_pal) +
# facet_wrap(~mod) +
facet_grid(rows = vars(mod), cols = vars(scenario),
switch = "y"))Code
plotly::ggplotly(ggplot(trace_tbl %>%
filter(state %in% c("BS_OAT_SS", "BS_SS")), aes(x = cycle_num, y = count, colour = state)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
ggrepel::geom_text_repel(data = trace_tbl %>%
filter(state %in% c("BS_OAT_SS", "BS_SS")) %>% filter(cycle_num == 180),
aes(x = cycle_num, y = count, label = state),
check_overlap = T, size = 3) +
scale_color_manual(values = colors_pal) +
# facet_wrap(~mod) +
facet_grid(rows = vars(mod), cols = vars(scenario),
switch = "y"))Code
plotly::ggplotly(ggplot(trace_tbl %>%
filter(state %in% c("BO_OD_ILLICIT", "BO_OD_RX",
"BO_MOD_BI", "BO_SEVERE_BI")), aes(x = cycle_num, y = count, colour = state)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
ggrepel::geom_text_repel(data = trace_tbl %>%
filter(state %in% c("BO_OD_ILLICIT", "BO_OD_RX",
"BO_MOD_BI", "BO_SEVERE_BI")) %>% filter(cycle_num == 180),
aes(x = cycle_num, y = count, label = state),
check_overlap = T, size = 3) +
scale_color_manual(values = colors_pal) +
# facet_wrap(~mod) +
facet_grid(rows = vars(mod), cols = vars(scenario),
switch = "y"))Code
plotly::ggplotly(ggplot(trace_tbl %>%
filter(state %in% c("BR_ILLICIT", "BR_OAT_INI",
"BR_OAT_MAINT", "BR_OAT_SS",
"BR_SS", "BR_OD_ILLICIT")), aes(x = cycle_num, y = count, colour = state)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
ggrepel::geom_text_repel(data = trace_tbl %>%
filter(state %in% c("BR_ILLICIT", "BR_OAT_INI",
"BR_OAT_MAINT", "BR_OAT_SS",
"BR_SS", "BR_OD_ILLICIT")) %>% filter(cycle_num == 180),
aes(x = cycle_num, y = count, label = state),
check_overlap = T, size = 3) +
scale_color_manual(values = colors_pal) +
# facet_wrap(~mod) +
facet_grid(rows = vars(mod), cols = vars(scenario),
switch = "y"))5 CBA results
5.0.1 Scenario 1
Code
cost_diff <- c(mod_nalx_1$total_net_present_cost - mod_basecase$total_net_present_cost,
mod_ss_1$total_net_present_cost - mod_basecase$total_net_present_cost,
mod_pres_guid_1$total_net_present_cost - mod_basecase$total_net_present_cost,
mod_all_int_1$total_net_present_cost - mod_basecase$total_net_present_cost)
cost_diff_per <- c(round((cost_diff/c(mod_basecase$total_net_present_cost))*100,2))
deaths_eff_tbl <- trace_tbl %>%
filter(scenario == "Scenario 1") %>%
filter(cycle_num == 180) %>%
select(-cycle_num) %>%
pivot_wider(names_from = mod, values_from = count) %>%
filter(state %in% c("BO_DEATH", "BO_OD_DEATH")) %>%
mutate(nalx_effect = Naloxone - `No Interventions`,
ss_effect = `Safer Supply` - `No Interventions`,
pg_effect = `Prescription Guidelines` - `No Interventions`,
all_effect = `All Interventions` - `No Interventions`,
nalx_effect_per = round(((Naloxone - `No Interventions`)/`No Interventions`)*100, 2),
ss_effect_per = round(((`Safer Supply` - `No Interventions`)/`No Interventions`)*100, 2),
pg_effect_per = round(((`Prescription Guidelines` - `No Interventions`)/`No Interventions`)*100, 2),
all_effect_per = round(((`All Interventions` - `No Interventions`)/`No Interventions`)*100, 2)) %>%
select(-c(Naloxone, `No Interventions`, `Safer Supply`, `Prescription Guidelines`, `All Interventions`))
total_death_oddeaths <- trace_tbl %>%
filter(scenario == "Scenario 1") %>%
filter(cycle_num == 180) %>%
select(-cycle_num) %>%
pivot_wider(names_from = mod, values_from = count) %>%
filter(state %in% c("BO_DEATH", "BO_OD_DEATH"))
od_deaths_diff <- c(deaths_eff_tbl$nalx_effect[1], deaths_eff_tbl$ss_effect[1],
deaths_eff_tbl$pg_effect[1], deaths_eff_tbl$all_effect[1])
od_deaths_diff_per <- c(deaths_eff_tbl$nalx_effect_per[1], deaths_eff_tbl$ss_effect_per[1],
deaths_eff_tbl$pg_effect_per[1], deaths_eff_tbl$all_effect_per[1])
deaths_diff <- c(deaths_eff_tbl$nalx_effect[2], deaths_eff_tbl$ss_effect[2],
deaths_eff_tbl$pg_effect[2], deaths_eff_tbl$all_effect[2])
deaths_diff_per <- c(deaths_eff_tbl$nalx_effect_per[2], deaths_eff_tbl$ss_effect_per[2],
deaths_eff_tbl$pg_effect_per[2], deaths_eff_tbl$all_effect_per[2])
saveRDS(cbind.data.frame(cost_diff, cost_diff_per,
od_deaths_diff, od_deaths_diff_per,
deaths_diff, deaths_diff_per), "point_estiamte.RDS")
saveRDS(total_death_oddeaths, "total_deaths_oddeaths.RDS")
costs <- paste0(round(cost_diff/1000000, 0), " (", cost_diff_per, "%)")
deaths <- paste0(round(deaths_diff, 0), " (", deaths_diff_per, "%)")
oddeaths <- paste0(round(od_deaths_diff, 0), " (", od_deaths_diff_per, "%)")
effects_tbl <- cbind.data.frame(mod_names[-1], costs, deaths, oddeaths)
effects_tbl |>
flextable() |>
add_header_row(colwidths = c(1, 3),
values = c("", "Mean Change Compared to No Intervention")) |>
set_header_labels(values = list(
`mod_names[-1]` = "Interventions",
costs = "Discounted Net Present Costs \n in Millions (%)",
deaths = "Deaths (%)",
oddeaths = "Overdose-related Deaths (%)")) |>
theme_vanilla() |>
set_table_properties(layout = "autofit") |>
align_text_col(align = "center", header = TRUE, footer = TRUE)Mean Change Compared to No Intervention | |||
|---|---|---|---|
Interventions | Discounted Net Present Costs | Deaths (%) | Overdose-related Deaths (%) |
Naloxone | 102 (0.01%) | 396 (0.01%) | -6420 (-4.05%) |
Safer Supply | 4233 (0.53%) | -13956 (-0.31%) | -3138 (-1.98%) |
Prescription Guidelines | -26103 (-3.28%) | 15252 (0.34%) | -1764 (-1.11%) |
All Interventions | -24876 (-3.13%) | 14115 (0.31%) | -8543 (-5.39%) |
5.0.2 Scenario 2
Code
cost_diff_2 <- c(mod_nalx_2$total_net_present_cost - mod_basecase_2$total_net_present_cost,
mod_ss_2$total_net_present_cost - mod_basecase_2$total_net_present_cost,
mod_pres_guid_2$total_net_present_cost - mod_basecase_2$total_net_present_cost,
mod_all_int_2$total_net_present_cost - mod_basecase_2$total_net_present_cost)
cost_diff_per_2 <- c(round((cost_diff_2/c(mod_basecase_2$total_net_present_cost))*100,2))
deaths_eff_tbl_2 <- trace_tbl %>%
filter(scenario == "Scenario 2") %>%
filter(cycle_num == 180) %>%
select(-cycle_num) %>%
pivot_wider(names_from = mod, values_from = count) %>%
filter(state %in% c("BO_DEATH", "BO_OD_DEATH")) %>%
mutate(nalx_effect = Naloxone - `No Interventions`,
ss_effect = `Safer Supply` - `No Interventions`,
pg_effect = `Prescription Guidelines` - `No Interventions`,
all_effect = `All Interventions` - `No Interventions`,
nalx_effect_per = round(((Naloxone - `No Interventions`)/`No Interventions`)*100, 2),
ss_effect_per = round(((`Safer Supply` - `No Interventions`)/`No Interventions`)*100, 2),
pg_effect_per = round(((`Prescription Guidelines` - `No Interventions`)/`No Interventions`)*100, 2),
all_effect_per = round(((`All Interventions` - `No Interventions`)/`No Interventions`)*100, 2)) %>%
select(-c(Naloxone, `No Interventions`, `Safer Supply`, `Prescription Guidelines`, `All Interventions`))
od_deaths_diff_2 <- c(deaths_eff_tbl_2$nalx_effect[1], deaths_eff_tbl_2$ss_effect[1],
deaths_eff_tbl_2$pg_effect[1], deaths_eff_tbl_2$all_effect[1])
od_deaths_diff_per_2 <- c(deaths_eff_tbl_2$nalx_effect_per[1], deaths_eff_tbl_2$ss_effect_per[1],
deaths_eff_tbl_2$pg_effect_per[1], deaths_eff_tbl_2$all_effect_per[1])
deaths_diff_2 <- c(deaths_eff_tbl_2$nalx_effect[2], deaths_eff_tbl_2$ss_effect[2],
deaths_eff_tbl_2$pg_effect[2], deaths_eff_tbl_2$all_effect[2])
deaths_diff_per_2 <- c(deaths_eff_tbl_2$nalx_effect_per[2], deaths_eff_tbl_2$ss_effect_per[2],
deaths_eff_tbl_2$pg_effect_per[2], deaths_eff_tbl_2$all_effect_per[2])
costs_2 <- paste0(round(cost_diff_2/1000000, 0), " (", cost_diff_per_2, "%)")
deaths_2 <- paste0(round(deaths_diff_2, 0), " (", deaths_diff_per_2, "%)")
oddeaths_2 <- paste0(round(od_deaths_diff_2, 0), " (", od_deaths_diff_per_2, "%)")
effects_tbl_2 <- cbind.data.frame(mod_names[-1], costs_2, deaths_2, oddeaths_2)
effects_tbl_2 |>
flextable() |>
add_header_row(colwidths = c(1, 3),
values = c("", "Mean Change Compared to No Intervention")) |>
set_header_labels(values = list(
`mod_names[-1]` = "Interventions",
costs_2 = "Discounted Net Present Costs \n in Millions (%)",
deaths_2 = "Deaths (%)",
oddeaths_2 = "Overdose-related Deaths (%)")) |>
theme_vanilla() |>
set_table_properties(layout = "autofit") |>
align_text_col(align = "center", header = TRUE, footer = TRUE)Mean Change Compared to No Intervention | |||
|---|---|---|---|
Interventions | Discounted Net Present Costs | Deaths (%) | Overdose-related Deaths (%) |
Naloxone | 85 (0.01%) | 303 (0.01%) | -4861 (-3.9%) |
Safer Supply | 4047 (0.51%) | -12743 (-0.28%) | -2435 (-1.95%) |
Prescription Guidelines | -26098 (-3.29%) | 15308 (0.34%) | -1485 (-1.19%) |
All Interventions | -24882 (-3.13%) | 14105 (0.31%) | -6607 (-5.3%) |
Code
tbl_effect_plot1 <- cbind.data.frame(Intervention = mod_names[-1], cost_diff_per,
deaths_diff_per, od_deaths_diff_per) %>%
pivot_longer(2:4, names_to = "group", values_to = "percent_change") %>%
mutate(scenario = "Scenario 1") %>%
bind_rows(., cbind.data.frame(Intervention = mod_names[-1], cost_diff_per_2,
deaths_diff_per_2, od_deaths_diff_per_2) %>%
pivot_longer(2:4, names_to = "group", values_to = "percent_change") %>%
mutate(scenario = "Scenario 2")) %>%
mutate(group = ifelse(group == "cost_diff_per_2", "cost_diff_per",
ifelse(group == "deaths_diff_per_2", "deaths_diff_per",
ifelse(group == "od_deaths_diff_per_2", "od_deaths_diff_per", group))))5.0.2.1 Scenario 1
Code
p7 <- ggplot(data = tbl_effect_plot1 %>%
filter(scenario == "Scenario 1"), aes(x = group, y = percent_change))+
geom_point(aes(colour = Intervention, shape = scenario), size = 1.5) +
scale_color_brewer(palette = "Dark2") +
scale_shape_manual(values = c(15, 16)) +
geom_hline(yintercept = 0) +
ylab("Percent Change Compared to No Intervention (%)") +
xlab("") +
scale_x_discrete(name ="",
labels=c("Costs","Deaths","OD-releated Deaths"))+
theme_minimal()
ggplotly(p7)5.0.2.2 Scenario 2
Code
p8 <- ggplot(data = tbl_effect_plot1 %>%
filter(scenario == "Scenario 2"), aes(x = group, y = percent_change))+
geom_point(aes(colour = Intervention, shape = scenario), size = 1.5) +
scale_color_brewer(palette = "Dark2") +
scale_shape_manual(values = c(15, 16)) +
geom_hline(yintercept = 0) +
ylab("Percent Change Compared to No Intervention (%)") +
xlab("") +
scale_x_discrete(name ="",
labels=c("Costs","Deaths","OD-releated Deaths"))+
theme_minimal()
ggplotly(p8)5.0.2.3 Both
Code
ggplotly(ggplot(data = tbl_effect_plot1, aes(x = group, y = percent_change))+
geom_point(aes(colour = Intervention, shape = scenario), size = 1.5) +
scale_color_brewer(palette = "Dark2") +
# scale_shape_manual(values = c(15, 16, 17, 18)) +
scale_shape_manual(values = c(15, 16)) +
geom_hline(yintercept = 0) +
ylab("Percent Change Compared to No Intervention (%)") +
xlab("") +
scale_x_discrete(name ="",
labels=c("Costs","Deaths","OD-releated Deaths"))+
theme_bw() +
facet_wrap(~scenario))Code
tbl_effect_plot2 <- cbind.data.frame(Intervention = mod_names[-1], cost_diff_per,
od_deaths_diff_per) %>%
mutate(scenario = "Scenario 1") %>%
bind_rows(., cbind.data.frame(Intervention = mod_names[-1], cost_diff_per_2,
od_deaths_diff_per_2) %>%
rename(cost_diff_per = "cost_diff_per_2",
od_deaths_diff_per = "od_deaths_diff_per_2") %>%
mutate(scenario = "Scenario 2"))5.0.2.4 Scenario 1
Code
p9 <- ggplot(data = tbl_effect_plot2 %>% filter(scenario == "Scenario 1"),
aes(x = cost_diff_per, y = od_deaths_diff_per))+
geom_point(aes(color = Intervention, shape = scenario)) +
scale_color_brewer(palette = "Dark2") +
scale_shape_manual(values = c(15, 16)) +
geom_hline(yintercept = 0) +
geom_vline(xintercept = 0) +
xlab("Mean Change in Costs Compared to No Intervention (%)") +
ylab("Mean Change in Overdose-related deaths \n Compared to No Intervention (%)") +
theme_minimal()
ggplotly(p9)5.0.2.5 Scenario 2
Code
p10 <- ggplot(data = tbl_effect_plot2 %>% filter(scenario == "Scenario 2"),
aes(x = cost_diff_per, y = od_deaths_diff_per))+
geom_point(aes(color = Intervention, shape = scenario)) +
scale_color_brewer(palette = "Dark2") +
scale_shape_manual(values = c(15, 16)) +
geom_hline(yintercept = 0) +
geom_vline(xintercept = 0) +
xlab("Mean Change in Costs Compared to No Intervention (%)") +
ylab("Mean Change in Overdose-related deaths \n Compared to No Intervention (%)") +
theme_minimal()
ggplotly(p10)5.0.2.6 Both
Code
ggplotly(ggplot(data = tbl_effect_plot2, aes(x = cost_diff_per, y = od_deaths_diff_per))+
geom_point(aes(color = Intervention, shape = scenario)) +
scale_color_brewer(palette = "Dark2") +
scale_shape_manual(values = c(15, 16)) +
geom_hline(yintercept = 0) +
geom_vline(xintercept = 0) +
xlab("Mean Change in Costs Compared to No Intervention (%)") +
ylab("Mean Change in Overdose-related deaths \n Compared to No Intervention (%)") +
theme_bw() +
facet_wrap(~scenario))